Intersecting raster and vector data
Clémentine Cottineau
Question: “How can I crop raster objects to vector objects, and extract the summary of raster pixels?”
Objectives:
- Crop a raster to the extent of a vector layer.
- Extract values from a raster that correspond to a vector file overlay.
Keypoints:
- Use the
crop()function to crop a raster object. - Use the
extract()function to extract pixels from a raster object that fall within a particular extent boundary. - Use the
extent()function to define an extent.
library(dplyr)
library(sf)
library(tibble)
library(ggplot2)
library(raster)
library(rgdal)
library(ggplot2)
library(dplyr)
library(here)Introduction
copied from the carpentry lesson Manipulating Raster Data).
We often work with spatial layers that have different spatial extents. The spatial extent of a shapefile or R spatial object represents the geographic “edge” or location that is the furthest north, south east and west. Thus is represents the overall geographic coverage of the spatial object.
The graphic below illustrates the extent of several of the spatial
layers that we have worked with in this workshop:
- Area of interest (AOI) – blue
- Roads and trails – purple
- Vegetation plot locations (marked with white dots)– black
- A canopy height model (CHM) in GeoTIFF format – green
Image Source: DCC
Frequent use cases of cropping a raster file include reducing file size and creating maps. Sometimes we have a raster file that is much larger than our study area or area of interest. It is often more efficient to crop the raster to the extent of our study area to reduce file sizes as we process our data. Cropping a raster can also be useful when creating pretty maps so that the raster layer matches the extent of the desired vector layers.
Import the raster
Data available here.
DSM_TUD <- raster(here("data","tud-dsm.tif"))
DTM_TUD <- raster(here("data","tud-dtm.tif"))
CHM_TUD <- DSM_TUD - DTM_TUD
CHM_TUD_df <- as.data.frame(CHM_TUD, xy = TRUE)
oai_boundary_tudlib <- st_as_sfc(st_bbox(raster(here("data","tudlib-rgb.tif"))))Crop a Raster Using Vector Extent
We can use the crop() function to crop a raster to the
extent of another spatial object. To do this, we need to specify the
raster to be cropped and the spatial object that will be used to crop
the raster. R will use the extent of the spatial object as the cropping
boundary.
To illustrate this, we will crop the Canopy Height Model (CHM) to
only include the area of interest (AOI). Let’s start by plotting the
full extent of the CHM data and overlay where the AOI falls within it.
The boundaries of the AOI will be colored blue, and we use
fill = NA to make the area transparent.
ggplot() +
geom_raster(data = CHM_TUD_df, aes(x = x, y = y, fill = layer)) +
scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) +
geom_sf(data = oai_boundary_tudlib, color = "blue", fill = NA) +
coord_sf()
Now that we have visualized the area of the CHM we want to subset, we
can perform the cropping operation. We are going to use the
crop() function from the raster package to
create a new object with only the portion of the CHM data that falls
within the boundaries of the AOI.
CHM_TUD_Cropped <- crop(x = CHM_TUD, y = st_as_sf(oai_boundary_tudlib))Now we can plot the cropped CHM data, along with a boundary box
showing the full CHM extent. However, remember, since this is raster
data, we need to convert to a data frame in order to plot using ggplot.
To get the boundary box from CHM, the st_bbox() will
extract the 4 corners of the rectangle that encompass all the features
contained in this object. The st_as_sfc() converts these 4
coordinates into a polygon that we can plot:
CHM_TUD_Cropped_df <- as.data.frame(CHM_TUD_Cropped, xy = TRUE)
ggplot() +
geom_sf(data = st_as_sfc(st_bbox(CHM_TUD)), fill = "green",
color = "green", alpha = .2) +
geom_raster(data = CHM_TUD_Cropped_df,
aes(x = x, y = y, fill = layer)) +
scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) +
coord_sf()
The plot above shows that the full CHM extent (plotted in green) is much larger than the resulting cropped raster. Our new cropped CHM now has the same extent as the aoi_boundary_HARV object that was used as a crop extent (blue border below).
ggplot() +
geom_raster(data = CHM_TUD_Cropped_df,
aes(x = x, y = y, fill = layer)) +
geom_sf(data = oai_boundary_tudlib, color = "blue", fill = NA) +
scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) +
coord_sf()
We can look at the extent of all of our other objects for this field site.
st_bbox(CHM_TUD) xmin ymin xmax ymax
83569.5 445251.5 87180.0 447180.0 st_bbox(CHM_TUD_Cropped) xmin ymin xmax ymax
85272.0 446295.0 85661.5 446694.0 st_bbox(oai_boundary_tudlib) xmin ymin xmax ymax
85272.00 446295.20 85661.28 446694.24 # plot_locations_HARV <-
# read.csv("data/NEON-DS-Site-Layout-Files/HARV/HARV_PlotLocations.csv")
# point_HARV <- st_read("data/NEON-DS-Site-Layout-Files/HARV/HARVtower_UTM18N.shp")
# utm18nCRS <- st_crs(point_HARV)
# plot_locations_sp_HARV <- st_as_sf(plot_locations_HARV, coords = c("easting", "northing"), crs = utm18nCRS)
#
# st_bbox(plot_locations_sp_HARV)
leisure_locations_selection <- st_read(here("data", "delft-leisure.shp")) %>%
filter(leisure %in% c("playground", "picnic_table"))Reading layer `delft-leisure' from data source
`/Users/ccottineau/Documents/GitHub/geospatial-data-carpentry-tud-2022-11/data/delft-leisure.shp'
using driver `ESRI Shapefile'
Simple feature collection with 298 features and 2 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 81863.21 ymin: 442621.1 xmax: 87370.15 ymax: 449345.1
Projected CRS: Amersfoort / RD Newst_bbox(leisure_locations_selection) xmin ymin xmax ymax
81863.21 442792.82 86719.87 449007.92 Our plot location extent is not the largest but is larger than the AOI Boundary. It would be nice to see our vegetation plot locations plotted on top of the Canopy Height Model information.
Challenge :
- Crop the Canopy Height Model to the extent of the study plot locations.
- Plot the vegetation plot location points on top of the Canopy Height Model.
Solution
CHM_plots_TUDcrop <- crop(x = CHM_TUD, y = leisure_locations_selection)
CHM_plots_TUDcrop_df <- as.data.frame(CHM_plots_TUDcrop, xy = TRUE)
ggplot() +
geom_raster(data = CHM_plots_TUDcrop_df, aes(x = x, y = y, fill = layer)) +
scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) +
geom_sf(data = leisure_locations_selection) +
coord_sf()
In the plot above, created in the challenge, all the vegetation plot locations (black dots) appear on the Canopy Height Model raster layer except for one. One is situated on the blank space to the left of the map. Why?
A modification of the first figure in this episode is below, showing
the relative extents of all the spatial objects. Notice that the extent
for our vegetation plot layer (black) extends further west than the
extent of our CHM raster (bright green). The crop()
function will make a raster extent smaller, it will not expand the
extent in areas where there are no data. Thus, the extent of our
vegetation plot layer will still extend further west than the extent of
our (cropped) raster data (dark green).
# Define an extent
So far, we have used a shapefile to crop the extent of a raster
dataset. Alternatively, we can also the extent() function
to define an extent to be used as a cropping boundary. This creates a
new object of class extent. Here we will provide the
extent() function our xmin, xmax, ymin, and ymax (in that
order).
# extent(CHM_TUD_Cropped_df)
new_extent <- extent(85272.25, 85661.25, 446295.2, 446693.8)
class(new_extent)[1] "Extent"
attr(,"package")
[1] "raster"TIP: The extent can be created from a numeric vector (as
shown above), a matrix, or a list. For more details see the
extent() function help file
(?raster::extent).
Once we have defined our new extent, we can use the crop() function to crop our raster to this extent object.
CHM_TUD_manual_cropped <- crop(x = CHM_TUD, y = new_extent)To plot this data using ggplot() we need to convert it
to a dataframe.
CHM_TUD_manual_cropped_df <- as.data.frame(CHM_TUD_manual_cropped, xy = TRUE)Now we can plot this cropped data. We will show the AOI boundary on the same plot for scale.
ggplot() +
geom_sf(data = oai_boundary_tudlib, color = "blue", fill = NA) +
geom_raster(data = CHM_TUD_manual_cropped_df,
aes(x = x, y = y, fill = layer)) +
scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) +
coord_sf()
Extract Raster Pixels Values Using Vector Polygons
Often we want to extract values from a raster layer for particular locations - for example, plot locations that we are sampling on the ground. We can extract all pixel values within 20m of our x,y point of interest. These can then be summarized into some value of interest (e.g. mean, maximum, total).
To
do this in R, we use the extract() function. The
extract() function requires:
The raster that we wish to extract values from, The vector layer
containing the polygons that we wish to use as a boundary or boundaries,
we can tell it to store the output values in a data frame using
df = TRUE. (This is optional, the default is to return a
list, NOT a data frame.) . We will begin by extracting all canopy height
pixel values located within our aoi_boundary_HARV polygon
which surrounds the tower located at the NEON Harvard Forest field
site.
tree_height <- extract(x = CHM_TUD, y = st_as_sf(oai_boundary_tudlib), df = TRUE)
str(tree_height)'data.frame': 621642 obs. of 2 variables:
$ ID : num 1 1 1 1 1 1 1 1 1 1 ...
$ layer: num 5.57 5.22 5.18 4.77 2.88 ...When we use the extract() function, R extracts the value
for each pixel located within the boundary of the polygon being used to
perform the extraction - in this case the aoi_boundary_HARV
object (a single polygon). Here, the function extracted values from
18,450 pixels.
We can create a histogram of tree height values within the boundary
to better understand the structure or height distribution of trees at
our site. We will use the column layer from our data frame
as our x values, as this column represents the tree heights for each
pixel.
ggplot() +
geom_histogram(data = tree_height, aes(x = layer)) +
ggtitle("Histogram of CHM Height Values (m)") +
xlab("Tree Height") +
ylab("Frequency of Pixels")
We can also use the summary() function to view
descriptive statistics including min, max, and mean height values. These
values help us better understand vegetation at our field site.
summary(tree_height$layer) Min. 1st Qu. Median Mean 3rd Qu. Max.
-3.460 0.000 0.375 4.343 8.523 36.729 Summarize Extracted Raster Values
We often want to extract summary values from a raster. We can tell R
the type of summary statistic we are interested in using the
fun = argument. Let’s extract a mean height value for our
AOI. Because we are extracting only a single number, we will not use the
df = TRUE argument.
mean_tree_building_height_AOI <- extract(x = CHM_TUD, y = st_as_sf(oai_boundary_tudlib), fun = mean)
mean_tree_building_height_AOI [,1]
[1,] 4.342554It appears that the mean height value, extracted from our LiDAR data derived canopy height model is 22.43 meters.
Extract Data using x,y Locations
We can also extract pixel values from a raster by defining a buffer
or area surrounding individual point locations using the
extract() function. To do this we define the summary
argument (fun = mean) and the buffer distance
(buffer = 20) which represents the radius of a circular
region around each point. By default, the units of the buffer are the
same units as the data’s CRS. All pixels that are touched by the buffer
region are included in the extract.
Image Source:National Ecological Observatory Network (NEON)
Let’s put this into practice by figuring out the mean tree height in
the 20m around the tower location (point_HARV). Because we
are extracting only a single number, we will not use the
df = TRUE argument.
point_Delft <- st_read(here("data", "delft-leisure.shp"))Reading layer `delft-leisure' from data source
`/Users/ccottineau/Documents/GitHub/geospatial-data-carpentry-tud-2022-11/data/delft-leisure.shp'
using driver `ESRI Shapefile'
Simple feature collection with 298 features and 2 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 81863.21 ymin: 442621.1 xmax: 87370.15 ymax: 449345.1
Projected CRS: Amersfoort / RD Newmean_tree_height_tower <- extract(x = CHM_TUD,
y = point_Delft,
buffer = 20,
fun = mean)
mean_tree_height_tower [1] NA 0.18382418 NA NA NA NA
[7] 0.14741693 NA NA 0.95534103 NA NA
[13] NA NA NA NA NA NA
[19] NA 1.08576819 1.77810457 1.50813700 NA NA
[25] NA NA NA NA 3.41709123 NA
[31] NA 5.44526759 NA NA NA NA
[37] NA NA NA NA NA NA
[43] NA NA NA NA 1.26241786 NA
[49] NA NA NA NA NA 3.63058446
[55] NA NA NA 2.64284311 NA NA
[61] NA NA 2.18177335 NA NA NA
[67] NA NA NA NA 1.71895883 NA
[73] 4.59002702 NA 9.22261513 3.71037319 3.67939371 2.67243673
[79] 4.01660438 NA NA NA NA NA
[85] NA NA NA 5.09733638 NA NA
[91] 3.60859259 NA NA NA 1.93923010 5.21877630
[97] NA NA NA NA NA NA
[103] NA NA NA 1.77163454 NA NA
[109] NA NA 1.70136037 2.58267291 4.21806988 NA
[115] NA NA NA NA NA NA
[121] NA NA 10.30380493 2.28110616 NA NA
[127] NA NA NA 3.64850905 NA 0.08129702
[133] NA 2.06486905 11.40954040 NA NA NA
[139] 2.08243861 1.11961589 6.40306065 6.41666083 6.71693856 5.15843022
[145] 4.27773571 NA NA NA NA NA
[151] NA NA NA NA NA NA
[157] NA 3.32952196 NA NA NA NA
[163] NA NA NA NA NA NA
[169] NA NA NA NA 2.07421889 NA
[175] 1.94442299 2.60783294 NA NA NA NA
[181] NA 7.78255215 NA NA NA 1.21956164
[187] 2.31018698 NA NA NA NA 7.29981885
[193] 2.72441063 NA NA NA 0.04019622 6.68842712
[199] 6.17018350 1.59578616 0.66848060 5.40900358 NA NA
[205] NA 2.25253693 NA NA 3.20768940 NA
[211] NA 0.25329162 NA NA NA 1.00615888
[217] 3.88431955 NA 10.92242381 NA NA NA
[223] NA NA 1.77282882 NA NA NA
[229] NA NA NA NA 1.34524428 1.75811156
[235] 1.87248210 1.37764249 1.50010280 2.11123471 2.40555998 1.10978271
[241] 0.85010793 2.65174185 2.63599304 2.40464816 NA NA
[247] NA NA NA 10.16022778 NA NA
[253] NA NA NA NA NA NA
[259] NA NA 0.61521977 NA 2.39372841 NA
[265] NA NA NA NA NA NA
[271] NA NA NA NA NA NA
[277] NA NA 5.96494518 2.29402613 NA NA
[283] 2.52794784 0.73166368 NA 1.58280219 4.89034569 0.07626293
[289] 0.40960883 3.79256709 NA NA NA 3.07805207
[295] NA NA NA NAChallenge: Extract Raster Height Values For Plot Locations
- Use the plot locations object (
plot_locations_sp_HARV) to extract an average tree height for the area within 20m of each vegetation plot location in the study area. Because there are multiple plot locations, there will be multiple averages returned, so thedf = TRUEargument should be used. Create a plot showing the mean tree height of each area.
Solution
leisure_locations_selection <- st_read(here("data", "delft-leisure.shp")) %>%
filter(leisure %in% c("playground", "picnic_table"))Reading layer `delft-leisure' from data source
`/Users/ccottineau/Documents/GitHub/geospatial-data-carpentry-tud-2022-11/data/delft-leisure.shp'
using driver `ESRI Shapefile'
Simple feature collection with 298 features and 2 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 81863.21 ymin: 442621.1 xmax: 87370.15 ymax: 449345.1
Projected CRS: Amersfoort / RD New# extract data at each plot location
mean_tree_height_plots_TUD <- extract(x = CHM_TUD,
y = leisure_locations_selection,
buffer = 20,
fun = mean,
df = TRUE)
# view data
mean_tree_height_plots_TUD ID layer
1 1 NA
2 2 0.95534103
3 3 NA
4 4 NA
5 5 NA
6 6 NA
7 7 NA
8 8 NA
9 9 NA
10 10 NA
11 11 1.08576819
12 12 1.77810457
13 13 NA
14 14 NA
15 15 NA
16 16 NA
17 17 NA
18 18 3.41709123
19 19 NA
20 20 NA
21 21 5.44526759
22 22 NA
23 23 NA
24 24 NA
25 25 NA
26 26 NA
27 27 NA
28 28 NA
29 29 NA
30 30 NA
31 31 NA
32 32 NA
33 33 NA
34 34 NA
35 35 NA
36 36 NA
37 37 3.63058446
38 38 NA
39 39 NA
40 40 NA
41 41 2.64284311
42 42 NA
43 43 NA
44 44 NA
45 45 2.18177335
46 46 NA
47 47 NA
48 48 NA
49 49 NA
50 50 NA
51 51 NA
52 52 1.71895883
53 53 3.67939371
54 54 NA
55 55 NA
56 56 NA
57 57 1.93923010
58 58 5.21877630
59 59 NA
60 60 NA
61 61 NA
62 62 NA
63 63 NA
64 64 1.77163454
65 65 NA
66 66 NA
67 67 NA
68 68 1.70136037
69 69 4.21806988
70 70 NA
71 71 NA
72 72 NA
73 73 NA
74 74 NA
75 75 NA
76 76 NA
77 77 NA
78 78 2.28110616
79 79 NA
80 80 NA
81 81 3.64850905
82 82 NA
83 83 0.08129702
84 84 NA
85 85 2.06486905
86 86 11.40954040
87 87 NA
88 88 NA
89 89 NA
90 90 2.08243861
91 91 1.11961589
92 92 6.40306065
93 93 6.41666083
94 94 6.71693856
95 95 NA
96 96 NA
97 97 NA
98 98 NA
99 99 NA
100 100 NA
101 101 NA
102 102 NA
103 103 3.32952196
104 104 NA
105 105 NA
106 106 NA
107 107 NA
108 108 NA
109 109 NA
110 110 NA
111 111 NA
112 112 NA
113 113 NA
114 114 NA
115 115 NA
116 116 NA
117 117 NA
118 118 NA
119 119 1.94442299
120 120 NA
121 121 NA
122 122 NA
123 123 7.78255215
124 124 NA
125 125 NA
126 126 1.21956164
127 127 NA
128 128 NA
129 129 7.29981885
130 130 2.72441063
131 131 NA
132 132 6.17018350
133 133 5.40900358
134 134 NA
135 135 NA
136 136 NA
137 137 NA
138 138 NA
139 139 NA
140 140 3.88431955
141 141 NA
142 142 NA
143 143 NA
144 144 NA
145 145 NA
146 146 NA
147 147 NA
148 148 NA
149 149 1.34524428
150 150 1.75811156
151 151 1.87248210
152 152 1.37764249
153 153 1.50010280
154 154 2.11123471
155 155 2.40555998
156 156 1.10978271
157 157 0.85010793
158 158 2.65174185
159 159 2.63599304
160 160 2.40464816
161 161 NA
162 162 NA
163 163 10.16022778
164 164 NA
165 165 NA
166 166 NA
167 167 NA
168 168 NA
169 169 NA
170 170 NA
171 171 NA
172 172 2.39372841
173 173 NA
174 174 NA
175 175 NA
176 176 NA
177 177 NA
178 178 NA
179 179 NA
180 180 NA
181 181 NA
182 182 NA
183 183 NA
184 184 NA
185 185 NA
186 186 NA
187 187 NA
188 188 5.96494518
189 189 2.29402613
190 190 NA
191 191 2.52794784
192 192 0.73166368
193 193 NA
194 194 3.07805207
195 195 NA
196 196 NA
197 197 NA
198 198 NA# plot data
ggplot(data = mean_tree_height_plots_TUD, aes(ID, layer)) +
geom_col() +
ggtitle("Mean Tree Height at each Plot") +
xlab("Plot ID") +
ylab("Tree Height (m)")
Summary and keypoints.
We have seen how to crop a raster to the extent of a vector layer and
how to extract values from a raster that correspond to a vector file
overlay. In short: - Use the crop() function to crop a
raster object. - Use the extract() function to extract
pixels from a raster object that fall within a particular extent
boundary. - Use the extent() function to define an
extent.